In [ ]:
import datetime

# Third-party
from astropy.io import ascii
import astropy.coordinates as coord
import astropy.units as u
import astropy.time as at
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
%matplotlib inline

In [ ]:
import astroplan
from astroplan import Observer, FixedTarget
from astropy.time import Time

In [ ]:
mdm = Observer.at_site("MDM", timezone="America/Phoenix")
t1 = Time(datetime.datetime(2015, 10, 14, 18, 0))
t2 = t1 + 12*u.hour
time_range = Time([t1, t2])

In [ ]:
def coords_in_rect(c, corner_c):
    if not c.frame.is_equivalent_frame(corner_c[0].frame):
        raise ValueError("Frame mismatch.")
    
    min_lon = corner_c[0].spherical.lon
    min_lat = corner_c[0].spherical.lat
    max_lon = corner_c[1].spherical.lon
    max_lat = corner_c[1].spherical.lat
    
    return ((c.spherical.lon > min_lon) & (c.spherical.lon < max_lon) & 
            (c.spherical.lat > min_lat) & (c.spherical.lat < max_lat))

How many GASS targets are there?


In [ ]:
css = ascii.read("/Users/adrian/projects/triand-rrlyrae/data/catalina.csv")
print(css.colnames)

In [ ]:
linear = ascii.read("/Users/adrian/projects/triand-rrlyrae/data/linear.csv")
print(linear.colnames)

In [ ]:
css_c = coord.SkyCoord(ra=css['ra']*u.deg, dec=css['dec']*u.deg, distance=css['helio_dist']*u.kpc)
lin_c = coord.SkyCoord(ra=linear['ra']*u.deg, dec=linear['dec']*u.deg, distance=linear['helio_dist']*u.kpc)

Filter by window on sky


In [ ]:
window_corners = [coord.SkyCoord(l=100*u.deg, b=-20*u.deg,frame='galactic'), 
                  coord.SkyCoord(l=220*u.deg, b=-10*u.deg,frame='galactic')]
css_ix = coords_in_rect(css_c.galactic, window_corners)
lin_ix = coords_in_rect(lin_c.galactic, window_corners)
print("{} CSS targets, {} LINEAR targets in this window.".format(css_ix.sum(), lin_ix.sum()))

In [ ]:
gass_tbl = css[css_ix]
gass = coord.SkyCoord(l=gass_tbl['l']*u.deg, b=gass_tbl['b']*u.deg, 
                      distance=gass_tbl['helio_dist']*u.kpc, frame='galactic')

Filter by distance


In [ ]:
ix = (gass.distance > 2.5*u.kpc) & (gass.distance < 8.*u.kpc)
gass = gass[ix]
gass_tbl = gass_tbl[ix]

In [ ]:
print("{} GASS targets".format(len(gass)))

In [ ]:
pl.figure(figsize=(10,6))
pl.scatter(gass.l.degree, gass.b.degree, c=gass_tbl['VmagAvg'], marker='o')
pl.xlim(220,100)
pl.ylim(-20,-10)

In [ ]:
fig,axes = pl.subplots(1,2,figsize=(12,5))
n,bins,pa = axes[0].hist(gass.distance, bins=np.linspace(2,10,15))
axes[0].set_xlabel("Radial dist. [kpc]")
axes[0].set_xlim(bins.min(), bins.max())

n,bins,pa = axes[1].hist(gass.galactocentric.represent_as(coord.CylindricalRepresentation).rho, bins=bins+8)
axes[1].set_xlabel("Cylindrical dist. [kpc]")
axes[1].set_xlim(bins.min(), bins.max())

axes[0].set_ylabel("Number RRL")

Red circle below is brightness limit of Catalina


In [ ]:
fig,ax = pl.subplots(1,1,figsize=(6,6),subplot_kw =dict(polar=True))

ax.add_artist(mpl.patches.Circle((-8.,0), radius=2.5, transform=ax.transData._b, facecolor='r', alpha=0.2))
ax.add_artist(mpl.patches.Circle((-8.,0), radius=0.5, transform=ax.transData._b, facecolor='y', alpha=1.))

gass_cyl = gass.galactocentric.represent_as(coord.CylindricalRepresentation)
css_cyl = css_c.galactocentric.represent_as(coord.CylindricalRepresentation)

ax.plot(css_cyl.phi.to(u.radian)[np.abs(css_cyl.z) < 5*u.kpc], 
        css_cyl.rho.to(u.kpc)[np.abs(css_cyl.z) < 5*u.kpc], 
        color='k', linestyle='none', marker='.', alpha=0.1)
# ax.plot(gass_cyl.phi.to(u.radian), gass_cyl.rho.to(u.kpc), 
#         color='k', linestyle='none', marker='o')
ax.scatter(gass_cyl.phi.to(u.radian), gass_cyl.rho.to(u.kpc), 
           c=gass_tbl['VmagAvg'], marker='o')

ax.set_rmax(20.0)
ax.grid(True)

ticks = [5,10,15]
ax.set_rticks(ticks)
ax.set_yticklabels(['{0:d} kpc'.format(x) for x in ticks])
ax.set_xlabel("Galactic Longitude", labelpad=15)
ax.tick_params(axis='y', labelsize=20)
# fig.savefig("/Users/adrian/papers/proposals/MDM-2015/GASS.pdf")

# ------
fig,ax = pl.subplots(1,1,figsize=(7,6))

gass_cyl = gass.galactocentric.represent_as(coord.CylindricalRepresentation)
css_cyl = css_c.galactocentric.represent_as(coord.CylindricalRepresentation)

ax.plot(css_cyl.rho.to(u.kpc)[np.abs(css_cyl.z) < 5*u.kpc], 
        css_cyl.z.to(u.kpc)[np.abs(css_cyl.z) < 5*u.kpc], 
        color='k', linestyle='none', marker='.', alpha=0.25)
cc = ax.scatter(gass_cyl.rho.to(u.kpc), 
                gass_cyl.z.to(u.kpc), 
                c=gass_tbl['VmagAvg'], marker='o')
ax.set_xlim(0,25)
ax.set_ylim(-12.5,12.5)
fig.colorbar(cc)

# ticks = [10,20]
# ax.set_rticks(ticks)
# ax.set_yticklabels(['{0:d} kpc'.format(x) for x in ticks])
# ax.set_xlabel("Galactic Longitude", labelpad=15)
# ax.tick_params(axis='y', labelsize=20)
# # fig.savefig("/Users/adrian/papers/proposals/MDM-2015/GASS.pdf")

In [ ]:
pl.hist(gass_tbl['VmagAvg'])

In [ ]:
gass_tbl['ID2015'] = ["GASS2015RR{0:d}".format(x+1) for x in np.arange(len(gass_tbl)).astype(int)]

In [ ]:
# ascii.write(gass_tbl, "/Users/adrian/projects/triand-rrlyrae/data/targets/gass.txt")

In [ ]:
# ascii.write(gass_tbl[['ID2015','ra','dec']], "/Users/adrian/projects/triand-rrlyrae/data/targets/gass_targets_2015_short.txt")#, format="ascii")

In [ ]:
gass_tbl_sex = gass_tbl.copy()

ra = coord.Longitude(gass_tbl_sex['ra']*u.deg)
gass_tbl_sex['ra_sex'] = ra.to_string(unit=u.hour, precision=5, sep=' ')

dec = coord.Latitude(gass_tbl_sex['dec']*u.deg)
gass_tbl_sex['dec_sex'] = dec.to_string(unit=u.degree, precision=5, sep=' ')

ascii.write(gass_tbl_sex[['ID2015','ra_sex','dec_sex']], "/Users/adrian/projects/triand-rrlyrae/data/targets/gass_targets_2015_short_sexa.txt")

Re-doing TriAnd selection


In [ ]:
window_corners = [coord.SkyCoord(l=100*u.deg, b=-35*u.deg,frame='galactic'), 
                  coord.SkyCoord(l=160*u.deg, b=-15*u.deg,frame='galactic')]
css_ix2 = coords_in_rect(css_c.galactic, window_corners)
print("{} CSS targets in this window.".format(css_ix2.sum()))

In [ ]:
triand_tbl = css[css_ix2]
triand = coord.SkyCoord(l=triand_tbl['l']*u.deg, b=triand_tbl['b']*u.deg, 
                        distance=triand_tbl['helio_dist']*u.kpc, frame='galactic')

Filter distance


In [ ]:
ix = (triand.distance > 15*u.kpc) & (triand.distance < 21*u.kpc)
triand = triand[ix]
triand_tbl = triand_tbl[ix]

In [ ]:
print("{} TriAnd targets".format(len(triand)))

In [ ]:
fig,axes = pl.subplots(1,2,figsize=(12,5))
axes[0].hist(triand_tbl['helio_dist'], bins=np.linspace(15,25,12))
axes[0].set_xlabel("Radial dist. [kpc]")
axes[0].set_xlim(11,25)

# axes[1].hist(gc_cyl_triand.rho, bins=8)
# axes[1].set_xlabel("Cylindrical dist. [kpc]")
# axes[1].set_xlim(17,24)

axes[0].set_ylabel("Number RRL")

In [ ]:
# fig,ax = pl.subplots(1,1,figsize=(8,8),subplot_kw =dict(polar=True))

# ax.add_artist(mpl.patches.Circle((-8.,0), radius=0.5, transform=ax.transData._b, facecolor='y', alpha=1.))

# ax.plot(all_gc_cyl.phi.to(u.radian)[sky_window], all_gc_cyl.rho.to(u.kpc)[sky_window], 
#         color='k', linestyle='none', marker='o', alpha=0.25)

# ax.plot(gc_cyl_triand.phi.to(u.radian), gc_cyl_triand.rho.to(u.kpc), 
#         color='k', linestyle='none', marker='o')
# ax.set_rmax(30.0)
# ax.grid(True)

In [ ]:
# triand2013_tbl = ascii.read("/Users/adrian/projects/triand-rrlyrae/data/TriAnd_RRL_26mar15.csv")
triand2013_tbl = ascii.read("/Users/adrian/projects/triand-rrlyrae/data/publication_data.csv")
triand2013 = coord.SkyCoord(ra=triand2013_tbl['ra']*u.deg,
                            dec=triand2013_tbl['dec']*u.deg)

In [ ]:
pl.hist(triand_tbl['VmagAvg'])

Match to old data


In [ ]:
idx, sep2d, _ = triand.match_to_catalog_sky(triand2013.galactic)

In [ ]:
match_ix = (sep2d < 5*u.arcsec)
match_ix.sum()

In [ ]:
triand2013_tbl

In [ ]:
names = list()
old_names = list()

for i,x in enumerate(sep2d):
    names.append("TriAnd2015RR{0:d}".format(i+1))
    if x < 5*u.arcsec:
        old_names.append(triand2013_tbl[idx[i]]['name'])
    else:
        old_names.append("--")
        

triand_tbl['ID2015'] = names
triand_tbl['ID2013'] = old_names
triand_tbl

In [ ]:
# ascii.write(triand_tbl, "/Users/adrian/projects/triand-rrlyrae/data/targets/triand1_targets_2015.txt")#, format="ascii")

In [ ]:
# ascii.write(triand_tbl[triand_tbl['ID2013']=='--'][['ID2015','ra','dec']], "/Users/adrian/projects/triand-rrlyrae/data/targets/triand1_targets_2015_short.txt")#, format="ascii")

In [ ]:
triand_tbl_sex = triand_tbl.copy()

ra = coord.Longitude(triand_tbl_sex['ra']*u.deg)
triand_tbl_sex['ra_sex'] = ra.to_string(unit=u.hour, precision=5, sep=' ')

dec = coord.Latitude(triand_tbl_sex['dec']*u.deg)
triand_tbl_sex['dec_sex'] = dec.to_string(unit=u.degree, precision=5, sep=' ')

ascii.write(triand_tbl_sex[triand_tbl['ID2013']=='--'][['ID2015','ra_sex','dec_sex']], "/Users/adrian/projects/triand-rrlyrae/data/targets/triand1_targets_2015_short_sexa.txt")

Brightness limit of Catalina

Claim is that it is V ~ 12.5


In [ ]:
from gary.observation import distance, apparent_magnitude
from gary.observation.rrlyrae import M_V

In [ ]:
distance(12.5 - M_V(-1.5)).to(u.kpc)

In [ ]:
apparent_magnitude(M_V(-1.5), 10.*u.kpc)

Exposure times


In [ ]:
from scipy.interpolate import InterpolatedUnivariateSpline

In [ ]:
pl.semilogy([14.5,16,17.5], [150,600,2400])
pl.xlim(14, 18)

In [ ]:
def Vmag_to_exptime(V):
    s = InterpolatedUnivariateSpline([14.5,16,17.5], np.log10([150,600,2400]), k=1)
    y = 10**s(V)
    return y

In [ ]:
triand_exptimes = [Vmag_to_exptime(V) for V in triand_tbl['VmagAvg']]
GASS_exptimes = [Vmag_to_exptime(V) for V in gass_tbl['VmagAvg']]

In [ ]:
len(triand_exptimes), len(GASS_exptimes)

In [ ]:
# print("TriAnd", (sum(triand_exptimes)*u.second).to(u.hour))
print("GASS", (sum(GASS_exptimes)*u.second).to(u.hour))

In [ ]: